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CL|' Abstract 

(D ■ 

'^ ' We propose a non-perturbative method for computing the renor- 

mahzation constants of generic composite operators. This method is 

^ . intended to reduce some systematic errors, which are present when 

H \ one tries to obtain physical predictions from the matrix elements of 

lattice operators. We also present the results of a calculation of the 
renormalization constants of several two-fermion operators, obtained, 
with our method, by numerical simulation of QCD, on a 16^ x 32 
lattice, at /3 = 6.0. The results of this simulation are encouraging, 
and further applications to four-fermion operators and to the heavy 
quark effective theory are proposed. 
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1 Introduction 

Renormalization of lattice operators is a necessary ingredient to obtain many 
physical results from numerical simulations, such as meson decay constants, 
form factors, structure functions, mixing amplitudes, etc. In this paper we 
study the renormalization of composite operators in lattice QCD. We pro- 
pose a method which avoids completely the use of lattice perturbation theory 
and allows a non-perturbative determination of the renormalization constants 
of any composite operator. For illustrative purposes, we consider some spe- 
cific applications to matrix elements of two-fermion operators. The approach 
is particularly useful in those cases, such as the scalar or pseudoscalar den- 
sities, where it is not possible to use the chiral Ward identities to determine 
the renormalization constants non-perturbativelyQ ||ll]-||^. Moreover, our pro- 
posal can be applied to many other cases, such as the renormalization of the 
four-fermion operators of the effective weak Hamiltonian, and to heavy-light 
currents in the Heavy Quark Effective Theory. Preliminary results have been 
presented in ref. [^. A similar attempt, limited to the divergent part of the 
AJ = 1/2 penguin-operator, can be found in ref. ||^. 

Renormalized lattice operators must correspond, in the limit of infinite 
lattice cut-off, to finite chirally covariant operators, which obey the same 
renormalization conditions as those in the continuum. With just a few ex- 
ceptions, lattice perturbation theory is used to evaluate the renormalization 
constants of lattice operators. The problem of mixing with lower dimensional 
operators requires however, a non-perturbative subtraction of all power di- 
vergences [|ll)|ll^[i]- Apart from this special case, the use of perturbation 
theory, in the computation of multiplicative constants and dimensionless 
mixing coefficients, is well justified, provided that the lattice spacing a is 
sufficiently small, i.e. a~^ ^ ^qcd- However, in those cases where the 
results from perturbation theory, obtained using the bare lattice coupling 
constant a^' = 5'o/4vr as an expansion parameter, can be checked using 
non-perturbative methods, there is generally a significant discrepancy. Sev- 
eral solutions to this problem have been proposed so far: 

• In order to improve the convergence of the perturbative series, the au- 



^In the following we will denote the method to determine the renormalization constants, 
using the Ward identities, as the Ward identity method. 



thors of ref. |T0| have proposed a "tadpole improved" perturbation 
theory. Nevertheless, ignorance of higher order contributions still rep- 
resents an important source of uncertainty in the extraction of physical 
results. 

• In some limited cases, corresponding to finite operators, namely the vec- 
tor and axial vector currents, and the ratio of pseudoscalar and scalar 
densities, a fully non-perturbative determination of the renormalization 
constants can be obtained with the use of chiral Ward identities |lI]-[0]. 

• One may fix non-perturbative renormalization conditions directly on 
hadronic matrix elements. This procedure was used in ref. 0, to 
subtract the divergent part of the A/ = 1/2 operator. The price is that 
one has to sacrifice a physical prediction, for any subtraction imposed 
on hadronic states. Thus, when several renormalization conditions are 
necessary, one looses much predictive power |l|, |^. 

Our proposal is to impose renormalization conditions non-perturbatively, 
directly on quark and gluon Green functions, in a fixed gauge, with given off- 
shell external states, with large virtualities. 

The method consists in mimicking what is usually done in perturbation 
theory. One fixes the renormalization conditions of a certain operator by 
imposing that suitable Green functions, computed between external off-shell 
quark and gluon states, in a fixed gauge, coincide with their tree level value. 
For example, if we consider the generic two-quark operator Or = "^/'r?/', we 
may impose the condition 

^r(p|Orb)|p2_^2 = (p|Orb)o, (1) 

where F is one of the Dirac matrices and {p\Oy\p)q is the tree level matrix 
element. The extension to more complicated cases, including four-fermion 
operators and operator mixing, is straightforward. This procedure defines 
the same renormalized operators, i.e. the same Wilson coefficient functions, 
in all regularization schemes, provided they are expressed in terms of the 
same renormalized coupling constant. However, the coefficient functions now 
depend on the external states and on the gauge, which must be specified. 



In principle this scheme completely solves the problem of large corrections 
in lattice perturbation theory, which are automatically included in the renor- 
malization constants. Matching of lattice operators to the corresponding 
ones, renormalized in a continuum renormalization scheme (for definiteness, 
the MS scheme), then requires continuum perturbation theory only. When 
next-to-leading corrections are known, the error on the physical matrix el- 
ements is of 0((a^ )^), in the continuum perturbative expansion, which is 
expected to have a better convergence. A continuum perturbative calcula- 
tion of the matching conditions is a common step in all approaches, standard 
or "tadpole improved" perturbation theory and non-perturbative renormal- 
ization. Our proposal can be applied to any composite operator, unlike the 
Ward identity method, which is limited to a few, albeit very important, 
cases^ 

The feasibility of using Green functions computed between quark states, 
in a fixed gauge, was already discussed in ref. 0, in the context of non- 
perturbative renormalization, using Ward identities. Indeed, in the case of 
the axial current, it was shown that the signal was rather good and the 
result was in agreement with the corresponding hadronic gauge invariant 
determination. 

As far as the renormalization conditions are concerned, our proposal is 
expected to work, whenever it is possible to fix the virtuality of the external 
states /i and to satisfy the condition Aqcd ^ /U ^ 1/a, in order to keep 
under control both non-perturbative and discretization effects. We stress 
that this requirement is common to all methods. The results of our numerical 
investigation are encouraging: they suggest that there does exist a "window" 
in n, where the method can be applied, at values of the lattice coupling 
constant used in current lattice simulations. 

The plan for the remainder of the paper is as follows. In the next section 
the basic formulae necessary to define our non-perturbative method are pre- 
sented, using as an example logarithmically divergent two-quark operators, 
such as the scalar and pseudoscalar densities. The discussion is generalised 
to all composite operators in section 3, and in the following section the im- 
plementation of the method in numerical simulations for two-quark operators 



•^It is still true that the statistical accuracy is in general slightly better with the Ward 
identities method, in those cases where it can be applied, see sec. o. 



is discussed. In section 5 we give some details about the perturbative evalua- 
tion of the renormalisation constants on lattices of finite volume. The values 
of the renormalisation constants obtained in a numerical simulation are pre- 
sented in section 6, and are compared to the results from perturbation theory. 
The final section contains our conclusions, a summary of the numerical re- 
sults and a discussion of possible further applications of the ideas presented 
in this paper. 



2 Non-perturbative renormalization conditions 

In this section, the non-perturbative renormalization scheme is presented in 
detail. We also discuss the conditions which must be satisfied for our method 
to be applicable. 

To facilitate the discussion, it is convenient to classify composite oper- 
ators into three main classes, according to their ultra-violet behaviour, as 

1. Finite operators: these include the vector and axial vector currents. 
Another example is the ratio of the scalar and pseudoscalar renormal- 
ization constants, Zs/Zp. For these cases, the Ward identity method 
is also applicable. 

2. Logarithmically divergent operators: this is a large family of op- 
erators, relevant to hadron phenomenology. It includes scalar and pseu- 
doscalar densities, four-fermion AF = 2 operators, some components 
of the energy-momentum tensor, etc. 

3. Power divergent operators: these operators are present when mix- 
ing with lower dimensional operators is possible. This happens in reg- 
ularizations which have an intrinsic mass scale, such as the lattice reg- 
ularization. Examples are often encountered in phenomenological ap- 
plications of lattice QCD, such as four-fermion operators relevant to 
A/ = 1/2 transitions or operators of order 1/m in the heavy quark 
effective theory (HQET) ]ll-|ll. 



For simplicity, we first consider the two-fermion operators of class 2., and 
extend the discussion to other operators in this class, and also to those in 



classes 1. and 3., in the next section. Thus, the formulae given below apply 
directly to the pseudoscalar and scalar densities. Throughout our discussion, 
we assume that discretization errors are negligible: in field theory language, 
this is equivalent to the statement that the renormalized Green functions do 
not differ appreciably from their values in the limit of an infinite ultra-violet 
cut-off. The discussion is presented in the limit of small quark masses and in 
Euchdean space-time. The discretisation of the quark action is assumed to 
be performed a la Wilson, characterized by explicit chiral symmetry breaking 
of the LagrangianQ. The extension to staggered fermions is straightforward. 

Let us consider the forward amputated Green function To (pa), of a two- 
fermion bare lattice operator 0{a), computed between off-shell quark states 
with four- momentum p, with p^ = /x^, and in a fixed gauge, for example the 
Landau gauge. Without gauge-fixing all Green functions computed between 
quark and gluon external states are zero. This is also true when the gauge is 
not completely fixed, e.g. in the Coulomb gauge. We define the renormalized 
operator 0{fi), by introducing the renormalization constant Zq 

0{fi) = Zo{fia,g{a))0{a). (2) 

Zq is found, by imposing the renormalization condition 

Zo(/ia,5'(a))^^^(/ia,5'(«))ro(pa)|p2=M2 = 1, (3) 

where Z^ is the field renormalization constant, to be defined below. This 
procedure defines a renormalized operator 0{fi) which is independent of the 



regularization scheme |T^-[T^. It depends, however, on the external states 
and on the gauge. This does not affect the final results, which, combined 
with the continuum calculation of the renormalization conditions, at any 
given order of perturbation theory, will be gauge invariant and independent 
of the external states. Let us specify the different quantities entering eq. (|^). 
To is defined in terms of the expectation valueQ of the non-amputated Green 
function Goipa), and of the quark propagator S{pa) 

To{pa) = —tT[Ao{pa)Po), (4) 



■^This includes the SW-Clover action Q. 

^ "Expectation value" means, as usual, that one averages the Green functions over the 
gauge field configurations, generated by Monte Carlo simulation, see sec. H. 



where 

Ao(pa) = Sipa)-'Goipa)Sipa)-\ (5) 

Po is a suitable projector on the tree-level operator: Pq = I {Po = 75) for 
the scalar (pseudoscalar) density. The factor 1/12 ensures the correct overall 
normalization of the trace (colourxspin=12). Projectors are very convenient 
when defining Green functions, particularly in the non-perturbative case. 
They have been extensively used in refs. fl^, |19| . Of course one can also use 
other definitions of Zq- 

1/2 
ZJ is the renormalization constant of the fermion field. It can be de- 
fined in different ways, some of which are equivalent perturbatively. Beyond 
perturbation theory, the most natural definition of Z^ is obtained from the 
amputated Green function of the conserved vector current V'-' . Indeed, one 
knows that for V'" the renormalization constant is equal to one: 

Zyl = Tyc X Z^^ = — tr (Ayc(pa)7/,) |p2=^2 x Z^^ = 1, (6) 

which implies 

Ztb = ^tr (Ayc(pa)7j |p2=^2. (7) 

Equations (0)^(0) completely define our method. In the remainder of 
this section, we discuss some important aspects concerning its applicability. 

In the above formulae, for simplicity, we have always considered only 
forward matrix elements. In general, one has the freedom to define the 
renormalization conditions at different external momenta p and p' {p 7^ p'). 
The virtualities of the quark states must be much larger than Aqcd- The 
reason is that, in order to obtain the physical result, we have to combine 
the matrix element of the renormalized operator 0(/x), with a Wilson coeffi- 
cient function. The latter is computed in continuum perturbation theory, by 

expanding in a^ at a scale of order /x. Thus, for the vahdity of this pertur- 
bative calculation, /i must be large. An important question in our program 
is whether it is possible to find, on the lattice, a scale /i which is sufficiently 
low, in order to have small 0{a) effects, and sufficiently large, in order to 
have small higher order corrections. The range of /i, where the above con- 
ditions are satisfied, depends on the value of Qq of the numerical simulation. 

6 



We stress that, if such a window does not exist, in current lattice simula- 
tions, an accurate matching of lattice operators to continuum ones becomes 
impossible, not only in our approach, but also with any other method. 

It might be expected that the condition /i ^ Aqcd ensures that perturba- 
tion theory is valid. When spontaneous symmetry breaking occurs however, 
as is the case in QCD, a large value of /i may not be enough, because of the 
presence of the Goldstone boson, the pion in our case. At low momentum 
transfers q = p — p', Green functions can receive a non-perturbative contribu- 
tion from the pion pole, which cannot be computed. The naive expectation 
is that this effect vanishes as 1/g^, at large q^. However, with fermion fields, 
this contribution is proportional to l/p"^ = l//i^, even when g^ = 0. The 
proof is given in the appendix. We conclude that a large value of /i solves 
the problem. At low values of /x^, however, for the pseudoscalar density and 
the axial current, we expect to find visible effects, see sees. | and ||. 

A danger to the non-perturbative scheme may come from the presence of 
Gribov copies. We do not address this problem here; it has been investigated 
in refs. pO| , |2T| . The results of ref. [^ indicate that, at least for the quan- 



tities of interest in the present study, the "Gribov uncertainty" , in current 
lattice simulations, is at most of the same order as the statistical error. 

One may imagine applying the renormalization conditions at very large 
values of /3, on lattices which have small physical volumes. This might seem 
to overcome the problem of the existence of the window in /x. For physical ap- 
plications however, it is necessary to evolve the renormalization constants to 
smaller values of f3. Non-perturbative contributions, which are not detected 
at large /3's, may then become important. This demonstrates the necessity, 
for any method, of the existence of a window Aqcd -C /x -C 1/a, where a is 
the lattice spacing used in numerical simulations of physical quantities. Only 
then are non-perturbative effects and lattice artefacts small simultaneously. 



3 Applications of the method to generic com- 
posite operators 

In the previous section, using the renormalization of the scalar and pseu- 
doscalar densities as a prototype, the general strategy of the non-perturbative 

7 



scheme has been presented. In the present section the necessary modifications 
needed in other cases are briefiy described. We will follow the classification 
of operators introduced in sec. 0. 



3.1 Finite operators 

In the case of vector and axial vector currents, we do not have the freedom of 
fixing the renormalization conditions in an arbitrary way: the renormaliza- 
tion conditions are acceptable only when they are compatible with the Ward 
identities. To be specific, one can apply eq. (H) to the local vector current 
in order to determine Zyh 

V _ _ ^^ (8) 



i yh 

We now show that this definition of Zy, which uses the general prescription 
given in the previous section, coincides with the one derived from the vector 
current Ward identity on quark states |jl], p2| p| 



% 



Zy.ky.{v^\v-\)\=-^[S-\v^\)-S-\v-h). (9) 






where AyL{p+q/2,p — q/2) is the amputated vertex with momentum transfer 
q. By applying d/dqp, at g = 0, on both sides of eq. (^, one gets 



Zy.{Ay.{p)+q, — Ay.ip+l,p-h\,=o)=-^J-S'\p). (10) 



The second term on the l.h.s. vanishes when g ^ 0. By tracing eq. (|10|) 
with 7p, we then recover the definition of ZyL of eq. (H). 

One is tempted to follow the same path for the local axial vector current, 
i.e. to define 

^A. = ^. (11) 

However this definition of Zj^l is not equivalent to the one obtained by the 
Ward identity 



?M 



9x , n-U Q^ 



Z^.A^.(p+ |,p- |)] = -^(^,S~\p+l) + S-\p-lH). (12) 
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For the remainder of this subsection we work in lattice units, setting a = 1. 



To see this we proceed as for the vector current, applying djdqp to both sides 
ofeq. (0). 



-K^=|;^"(^)-|;^"M- (13) 

In the axial case, however, the second term on the l.h.s. does not vanish, as 
g — i> and, indeed, this term is necessary to saturate the Ward identity, in 
presence of a massless Goldstone boson. As demonstrated in the appendix 
however, this term becomes negligible for high values oi p^. In this case, by 
tracing eq. (]I3|), one recovers (|11|). In conclusion, we expect to find a value 
of Zj^L close to those obtained from the Ward identity on hadron or on quark 
states 0, by using eq. ([Tl|), at large /i^. 

Zs/Zp is another finite quantity, which is completely fixed by the Ward 
identities, even though the renormalization conditions on Zs and Zp are 
separately arbitrary. In other words one is free to decide the renormalization 
conditions of one of the two operators, the scalar density say, and the value 
of Zp will automatically follow. The procedure of sec. ^ respects the Ward 
identities and thus guarantees that the correct value of Zs/Zp is obtained. 
This is clearly true only at large values of /x^, for the same reasons as for the 
axial vector current. 



3.2 Logarithmically divergent operators 

For this category of operator only the question of operator-mixing is left to 
be discussed. There are two kinds of mixing: the one which can occur also in 
the continuum and the one induced by the explicit lattice chiral symmetry 
breaking. We are only concerned here with the latter case, and hence restrict 
our discussion to operators which renormalize multiplicatively in the contin- 
uum. To illustrate our arguments, we consider the four-fermion operator 



which is relevant to the calculation of the K^-K^ mixing amplitude. Because 
of the Wilson term, this operator mixes with operators of different chiralities 



Otif = Zo{O^J^' + Y.Z^O^)^ (15) 

i 

where the operators O* are given by O^ = {sd){sd),0^ = (s75(i)(s75(i), etc. 
2^ . The mixing coefficients Zi are finite functions of glia). The over- all 



renormalization Zq is necessary to eliminate the logarithmic divergence. 

We now sketch the renormalization procedure for O^^^^. Let us start 
from the four-point amputated, Green function of 0^^^^^, between quark 
states with momentum p^ = /i^ 

Af,T(H, (16) 

where a, f3 . . . and A,B... are spin and colour indices respectively. To find 
the mixing coefficients Z^, one uses suitable projectors P, as was done for 
two-quark operators^, for example 

tr(AP) = Af,^5^(75)/3„(75)... (17) 

By projecting on all the possible colour-spin structures, corresponding to 
the operators O*, one obtains a system of linear equations in the mixing 
coefficients Z^. Using the Z^, determined in this way, one can compute the 
Green function of the subtracted operator O^ = 0(^^^^^ + J2i ZiO^. From the 
amputated vertex of O^, one can finally remove the logarithmic divergence, 
by imposing the renormalization condition 

Zo(/i,^(a))z^^(/i,^(a))r^(pa)|p2=^2 = 1, (18) 

where 

r^~tr(A''Pz.xL), (19) 

and Plxl is the suitable projector on the 7^(1 — 75) ®7''(1 — 75) Dirac tensor. 

^ In order to simplify the presentation, we do not discuss here the comphcations arising 
from the different possible colour structures. 
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Figure 1: Diagram which mixes the four fermion operator in eq. pO| ) with 
the scalar and pseudoscalar densities. The mixing is induced by the chiral 
breaking term present in the Wilson and Clover actions. 

3.3 Power divergent operators 

Composite operators are divergent in 1/a when they can mix with lower 
dimensional operators. Examples are given by the operators relevant to deep 
inelastic scattering and by the penguin-operators appearing in the AS = 1/2 
Hamiltonian. Let us consider the AS = 1/2 operator 

OA5=l/2 ^ (-^^^1 _ ^^)^) J:i9^ril + 75)g.). (20) 

i 

Oy^^^~ mixes with the scalar and pseudoscalar denisties, even at zeroth 
order in ag, through the diagram in fig. [iQ 

OtT^' = ZoOtt'^' + Z,{-sd) + Z,{-s^,d) + .... (21) 

To find the coefficients 2'i2, one imposes that the two-point Green func- 
tion of the operator on the r.h.s. of eq. (^Tj), on off-shell quark states, is 
zero 0. Once all the divergent parts have been removed, one can impose 
the renormalization conditions, necessary to remove the over-all logarithmic 
divergence, as discussed above. 



^ '-^latt also mixes with the chromo-magnetic operator, e.g. scr^^G^'^d, which is not 
being considered here. 
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4 Implementation of the method in numeri- 
cal simulations 

To implement the ideas of the previous section, let us consider the generic 
local two-fermion operator 

Or(x) = ip{x)Ti;{x) (22) 

where F is any Dirac matrix. The lattice, bare Green function of Or, between 
off-shell quark states, can be obtained from 

1 ^ 
Go{x, y) = {ij{x)Orm{y)) = - J] S,{x\0)TS,{0\y), (23) 

^^ 1=1 

where i = 1,N labels the gauge field configurations and Si{x\0) is the quark 
propagator on the single configuration i, obtained by inverting the discretized 
Dirac operator. We denote by Si{x\y) the inverse Dirac operator, which is not 
translationally invariant, in contrast to S{x — y), which is the usual quark 
propagator. The gauge field configurations are generated with some stan- 
dard, gauge-invariant algorithm; the gauge fields and the quark propagators 
are then rotated into the Landau gauge, which will be defined in sec. ^ 
From eq. (|2|), using Si{x\y) = 755'i(?/|a;)^75, one gets 

Goipa) = f d^xd^ye-'^P-''-Py^Go{x,y) 

= ]vE(y d''xS,{x\0)e-^^-^)T{J dV,(0||/)e^^-^) 
1 ^ r 

1 ^ 
Y.S^{pm{J,Sj{p\Oh,), (24) 



Nt, 



where we have defined 

S^{p\y)= [ d''xS,{x\y)e-'^^--\ (25) 
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because, on a single configuration, the quark propagator is not translation- 
ally invariant. Traslational invariance is recovered after averaging over the 
configurations. We can then write 



1 ^ 
S{pa) = {S{pm = -J:S^{p\Q). (26) 

Eqs. ( plD and ( ^6]) give Goipo) and S{j)a) in terms of quark propagators, 
computed by inverting the Dirac operator in the presence of a background 
gauge field configuration, generated by numerical simulation. With these 
quantities at hand, one can apply the strategy and the formulae, introduced 
in sec. |^. We add some information, that can help the reader in practical 
applications. The quark wave-function renormalization constant Z^ has been 
computed using 

Z^ = TyL X ZyL = — tr (A^L (^0)7^) |p2=^2 X ZyL (27) 

where AyL^pa) is defined as in eq. (H); the trace is taken over colours and 
spins and ZyL is the renormalization constant of the local vector current 
V^(a;) = ip{x)'yfj,4'{x). Equations (|^) and (^71) are two equivalent definitions 
of Z^. We have used eq. ( P7D instead of eq. (J^), because the former is more 
convenient when, as in our case, one has only computed quark propagators 
which stem from the origin. The three-point {q — V — q) Green function 
is obtained by inserting the local vector current in the origin. ZyL can 
be obtained with very high accuracy from ratios of matrix elements of the 
conserved and local currents P, E3] . From the amputated Green functions 



ri(pa) = — tr(Ai(pa)i 
r75(P«) = T^tr(A^3(pa)75) 



r7^75(P«) = ^tr(A^^^3(pa)757^), (28) 

by imposing the renormalization conditions (^, at several values of p^ = yU^, 
to be specified in sec. P, we have determined Z^p^^f/x, (^(a) j, as functions of 



/.^ 
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If we wish to compute Zyh using the procedure defined in sec. || we 
clearly cannot take Z^ from eq. (^). From the Ward identity we have 

Z^ = -^_tr(— ^)|,.=,. (29) 

where dS{p)~^ j dj) = 1/4 x '~ipdS{p)~^ /dpp. To avoid derivatives with respect 
to a discrete variable, we have used 



/ X tr(-zEA=i.4 7Asin(pAa)g \pa)) 

In one loop perturbation theory, ZL = Z^, in the Landau gaugeQ. By using 

r7MM = ^tr(A^^(pa)7M), (31) 

and imposing the renormalization conditions (Q), we have determined ZyL at 
several /i^. From the general discussion of the previous section, one expects 
that ZyL, defined through the present procedure, is independent of /i^, up to 
terms of 0{a), as has been effectively found, see sec. ||. 

In actual numerical simulations, the renormalization procedure could be 
obscured by lattice artefacts, i.e. terms of 0{fia), which are present when 
one imposes the renormalization conditions at large values of /i^, required to 
avoid large non-perturbative or higher order effects. As shown in ref . [^ (we 
use here the same notation), in order to reduce discretization errors, due to 
the finite value of the lattice spacing, one can compute correlation functions 
using a nearest-neighbour improved fermion action [|T^] 

where Sp' is the usual Wilson action, r is the Wilson parameter and F^i, is the 
lattice field strength tensor. Moreover, it is also necessary to use improved 
operators Of defined as 

0^{x)=tlj{x)'R,Tl,tlj{x), (33) 



^In general, they differ by a finite term of order ag 
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where the rotations of the fermion fields are defined as 

-n.= {l-J-0^-m,)) , X={l + a-^(^y + mo)). (34) 

By using the improved action and operators, one ehminates discretization 
errors of 0{a), and only those of 0{asa) or 0{a^), or higher, remain. In eq. 
(|5HD, Kx and Ry can be replaced by 

^^ = l-J-'^^ , ^, = l + a^5^, (35) 

in computations of the hadron spectrum and of on-shell operator matrix ele- 



ments. As discussed in refs. [|26|, ^ , the improvement procedure is equivalent 
to expressing all correlation functions in terms of "effective" quark propaga- 
tors 

Sf\x\Q) = (l - a-'^.^Sl{x\Q)[l + a^5,=o) (36) 

where Sl{x\Q) is the fermion propagator of the improved theory from the 
point X to the origin. In order to eliminate terms of 0{a) in the off-shell 
correlation functions, used in this study, a further step is necessary 0: we 
have to shift the quark propagator by a local term, which in Fourier space 
appears as a constant, to be added to its diagonal components 

In the following, the vector current V^ , the axial vector current A^, the 
pseudoscalar density P and the scalar density 5* are improved operators, 
obtained by using eq. (^6|). 



5 Lattice perturbation theory 

In order to gain some insight into the non-perturbative results presented in 
the next section, we have performed the corresponding one-loop perturbative 
calculation both in the continuum and on the lattice in the Landau gauge. 



This calculation is an extension of the one presented in ref. |27 . 
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The standard calculation of the one-loop corrections corresponds to a 
computation of the renormalization constants with a fixed ultra-violet cut- 
off on the loop momenta, \pi\ < it /a, neglecting terms of 0{a) in the final 
result. The renormalization constants obtained in this limit depend on the 
quark mass rrig and on the external momenta only through the logarithmi- 
cally divergent terms, which are related to the anomalous dimension of the 
operator. Since, at one loop, these terms are universal, i.e. they are the same 
on the lattice and in any continuum regularization, the renormalization con- 
stants, necessary to relate lattice to continuum renormalized operators, are 
independent of the quark mass and external momenta (they depend however 
on the lattice spacing and on the continuum renormalization scale /i). We call 
this limit "Standard Perturbation Theory" (STP). The non-perturbative cal- 
culation however, is performed on a lattice of finite size, and at a fixed value 
of the lattice spacing and quark mass: thus we cannot make terms of 0{a) 
arbitrarily small. In particular, since we want to study the /^-dependence of 
the renormalization constants, it is important to know at which values of /x 
(and ruq), the perturbative results, obtained on a finite lattice and with a 
fixed lattice spacing, differ from those obtained in SPT. In order to follow 
closely the non-perturbative computation, the perturbative results have also 
been obtained, on a lattice of size 16^ x 32. In this way, we keep all terms 
of 0{asa) or O(a^), or higher, in the final result (we call this procedure 
"Discrete Perturbation Theory", DPT). The loop integrals were computed 
numerically, both in SPT and in DPT. 

The non-perturbative results have been obtained by tracing suitable am- 
putated Green functions, computed between off-shell quark states, in the 
Landau gauge, cf. eq.(|). In perturbation theory, this procedure differs from 
the standard one, where one identifies the correction from the term propor- 
tional to the tree-level matrix element of the operator under consideration 



Ig, |19|. The difference is due to finite terms which, in the standard proce- 
dure, are considered as part of the matrix element of the operator, whereas 
with the projectors, are considered as part of the one-loop corrections. At 
one-loop order, these terms are the same in the continuum and on the lattice 
(in the limit a — > 0) and cancel when one computes the difference between 
the continuum and the lattice renormalization constants. To be more spe- 
cific, let us consider the diagrams in fig. ^ computed in Naive Dimensional 
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a) 





Figure 2: Vertex diagrams with the Clover action. The second and third 
diagrams refer to the contribution from the rotated part of the operator. 




Figure 3: Self energy diagram. 



Regularization (NDR), in the Feynman gauge 



P ^ 

A7^(-) =< p\1p7^^'^P\p >-- 



^-S-K^ 






(37) 



where Cp = {N'^ — 1)/2A^, 1/e = 1/e — ln47r — je and p is the external 
momentum. Usually, one takes the term proportional to (1/e + 1), i.e. the 
coefficient of 7^, as the one loop contribution to the bare operator. With the 
projector (see eq.(^)) the last term in eq. (|37|) also gives a contribution, so 
that the factor now becomes (1/e + 1/2). In table |l|, we give the one-loop 
corrections obtained, with the projectors and in the 't Hooft and Veltman 
(HV) dimensional renormalization scheme |^, for the two-quark operators 
considered in this paper. The NDR regularization differs from the HV one 
by the definition of the ant i- commuting properties of 75. For details see for 
example [^ . In NDR the corrections to the axial current and pseudoscalar 
density coincide with those to the vector and scalar density respectively. 
By writing the gluon propagator as G{k) = {—d^yk"^ + (1 ~ ^)k^ky) //c^, we 
denote as "Landau" the contribution to the one-loop corrections proportional 
to (1-A). 

With the MS procedure, after the subtraction of the pole, the irreducible 
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Operator 


Feynman 


Landau 


7/. 


1/e +1/2 


-1/e -1/2 


7m75 


1/e + 5/2 


-1/e -1/2 


1 


4/e + 6 


-1/e - 2 


75 


4/e + 10 


-1/e -2 


S 


-1/e - 1 


1/e + l 



Table 1: The results of the calculation of the one-loop diagram of fig. ^for 
different two-quark operators, in the HV scheme, are given. We also give the 
result of the calculation of the self-energy diagram of fig. |^. The mass of the 
quark has been taken to be zero, the inverse quark propagator is written as 
ifL{p^) and a factor {ag/ 4:T{)C pijP' / n^j is omitted. 



vertex of fig. ^ has the formj^ 



To{p/^^) = 1 + f^C7p(-7(^) HpVf^') + c^')- 



(38) 



From table |I], one finds C^^ = 1/2, C^^ = 6, etc., in the Feynman gauge 
-1/2, C^^ = —2, etc., for the Landau (longitudinal) terms. 
On the lattice, one gets a similar expression 



and C^s 



Toipa) = 1 + f^C^(-7(^) Hp'a') + C^^ 



Latt\ 



(39) 



The values of Cjf"^^, for the Feynman and Landau corrections are given in 
table 0. In table ^ we have denoted by "Rotated Feynman" the contribution 
coming from the rotated part of the two-fermion operators, see eqs. (|33|) and 
([341). The Landau contribution to the rotated part of the operators vanishes. 
This is required by the fact that the renormalization constants which relate 
lattice to continuum operators are gauge invariant. 

Following refs. p^-p^, we write 



0(/i) 



Oo 



An 



-C^(-7(°)ln(/ia)2 + Ao) 0( 



(40) 



^To(pI p^) is obtained from Ko{p/ 1^) using a projector, as described in eq. 



Operator 


Feynman 


Rotated Feynman 


Landau 


7m 


6.62 


-3.52 


-4.29 


7m75 


5.09 


-11.68 


-4.29 


1 


16.11 


-12.00 


-5.79 


75 


19.18 


4.31 


-5.79 


S 


8.21 





4.79 



Table 2: The values of Cjf'^^ for different two-quark operators and for the 
self- energy (C^), with the Clover action. The mass of the quark has been 
taken to be zero. 



where Aq can be computed from tables |l] and | 






/^Latt 



/^MS 



y^Latt 



(41) 



In the difference between continuum and lattice one-loop corrections, all 
Landau pieces cancel, as expected. The results for A^^ , A^^^^ , ■ ■ ■ all agree 
with those of ref. ^7 . 



On a lattice, for a given quark mass and external momenta, there are 
terms of order a which can become important at large values of m^ and p. 
Since with our method we have to impose the renormalization conditions at 
values of the momenta which are large, it is important to monitor, at least in 
perturbation theory, the values of masses and momenta at which 0{a) effects 
become important. The Green function in DPT has a form similar to (|39D 



Toipa) 



'-> 



(-7(0) ln(/a2) + C^-'\pa, m,a)) 



(42) 



where, however, the "constant" Cq^** now depends on the terms of 0{a). 
Thus, in the calculation of the one-loop corrections in DPT, all the cut-off 
dependent terms have been included. On a finite volume, the propagators 
of massless particles at zero momentum must be regulated. The results 
of Feynman diagrams depend on the regulator, however this dependence 
vanishes as the volume is increased. In the results presented below, we have 
set the zero-momentum terms in both the quark and gluon propagators to 
zero. We observe a small difference in the results from SPT and DPT, at low 
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values of /x^, which is due to this choice. At the value of the quark mass at 
which we have performed the simulation, mass corrections are negligible. 

The coupling constant appearing in eqs. (^) and (^) is the bare lattice 
coupling. In ref. |10| it was argued that it is possible to optimize the con- 



vergence of the perturbative series by reorganizing the expansion in terms 
of a "boosted" coupling a^ , defined in terms of a physical quantity. In the 
comparison with the non-perturbative results below we follow the sugges- 



tion of ref. [|I^, and compute the perturbative corrections using an effective 
coupling 

«r = nr^lr , c.f^^^ - 1.68 af^^^, (43) 

where af^"^"^ = 5'o/4vr and {^Trllpiaq) is the expectation value of the pla- 
quette. We call this expansion "Boosted Discrete Perturbation Theory" 
(BDPT). Analogously we denote the "Boosted Standard Perturbation The- 
ory" by BSPT. 

The results of the calculation of the renormalization constants in BDPT 
will be compared to the non-perturbative results in the next section. 



6 Non-perturbative results 

In this section, we give the main results of our numerical study of the non- 
perturbative renormalization procedure proposed in this paper and a com- 
parison of these results with perturbation theory. 

We have performed a simulation, by generating 36 independent gluon 
field configurations, on a 16^ x 32 lattice, at /3 = 6.0. The quark propagators 
have been computed for a single value of the quark mass {arriq ~ 0.07), 
corresponding to a value of the hopping parameter K = 0.1425. All the 
Green functions have been computed in the lattice Landau gauge, which is 
obtained by minimizing the functional 

tT[Y^{u,ix) + Ulix))]. (44) 



Possible effects from Gribov copies or spurious solutions 0] have not been 



considered. As shown below, for those quantities which can be determined 
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also in a gauge invariant way, the non-perturbative results obtained on quark 
state in a fixed gauge are in good agreement with those obtained using the 
Ward identity method 0. 

Vector current: In fig. |, the renormalization constant of the local 
vector current ZyL, obtained using eqs. ( ^0|) and (|3T1) , is given. As expected, 



ZyL is independent of the scale, whithin statistical errors, up to large values 
of /i^, where distortions due to lattice discretization become important. This 
is a consequence of the equivalence, which can be established up to terms of 
O(asa), of the method used here to determine ZyL, and of the Ward identity 
for the local vector current, see sec. ^. By using the points at jji^a? ~ 1, 
corresponding to /i ~ 2 GeV, we get Zy^ = 0.84(1), to be compared with 
Z^l = 0.824(2) from the Ward identity method, and Z^^ = 0.83 (0.86-0.87 
on a 16^ X 32 lattice, at /x^a^ ~ 1) from BSPT. The three methods are in 
good agreement. 

Scalar density.: Zs is an ideal quantity with which to check the validity 
of our method for a number of reasons: it cannot be determined using the 
Ward identity method, it is logarithmically divergent and gauge dependent, 
and it is not affected by the pole of the Goldstone boson (in contrast to the 
pseudoscalar density and the axial current), sec. 0. 

In fig. ^, Zs is shown as a function of the renormalization scale. We 
report separately the results, obtained by using Z^, eq.(^7p, and Z'^, eq.(|5D|). 
Notice that the points, corresponding to different definitions of Z^, are in 
very good agreement where /i^a^ is not too large, and discretization errors 
are small. This can also be seen by comparing the results in BSPT (dashed 
curve) to those in BDPT (continous curve). We do not expect agreement with 
perturbation theory, either at low /x^, where higher order effects become very 
large, or at high /i^, where lattice distortions are important, as shown by the 
continous curve. The numerical results follow the theoretical expectations 
and we find good agreement between the non-perturbative determination of 
Zs and the predictions from boosted perturbation theory for 0.3 < n'^a'^ < 1 
(1.1 GeV < /i < 2 GeV). Notice the Zs (Zp) is defined here from the off- 
shell Green function and depends on the anomalous dimension of the scalar 
(pseudoscalar) density, and hence on log(/i^a^). Moreover, the result is gauge 
dependent, even in the continuum, as can be shown by an explicit calculation 
in one-loop perturbation theory. 
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Axial vector current: The local axial- vector current shares some of 
the features of the vector current. Its renormalization constant is finite at all 
orders in perturbation theory ^1\ and can be determined from Ward identi- 



ties. However, the axial-vector current is coupled to the would-be Goldstone 
boson of QCD, which can give an important contribution at low /x^. Zj^l, 
which like ZyL, should be independent of /i, is shown as a function of /i^ in 
fig-i 

We interpret the strong /i-dependence of Z^^l, at low /i^, as the non- 
perturbative effect of the pseudoscalar state. Unfortunately, there is no clear 
sign of the existence of a plateau between the non-perturbative regime and 
the large /i region, where lattice artefacts become important|^. Without the 
information coming from the Ward identity method |^, it would be difficult 
to determine Z^l confidently. Nevertheless, we do observe that, at fi'^a^ ~ 1, 
just before lattice artefacts become large in DPT, the results are close to the 
value determined with the Ward identity, Z^^ = 1.06(2) 0, ^, see fig. p. As 



observed in ref. P], perturbation theory gives values of Zj^l smaller than one, 
while the Ward identity method gives values larger than one, for f3 = 6.0-6.2. 
A value larger than one is also suggested by our non-perturbative results. 

Pseudoscalar density: The pseudoscalar density shares the main fea- 
tures of the scalar density, with two important differences. Firstly it is 
coupled to the would-be Goldstone boson, and secondly the one-loop per- 
turbative corrections, with the Clover action, are quite large. It is then not 
surprising that the non-perturbative value of the corresponding renormal- 
ization constant Zp lie well below the perturbative result, as shown in fig. 
1^. Had we used standard perturbation theory in fig. |^, instead of boosted 
perturbation theory, the discrepancy would have been even worse. The dis- 
crepancy could have been anticipated from the results of ref. 0, where it 
was shown that the ratio Zs/Zp, determined by using the Ward identity 
method, was significantly larger than the result obtained in boosted pertur- 
bation theory. Combining this information with the agreement between the 
non-perturbative determination of Zs and the result from BPT, one would 
conclude that the difference is due to Zp. This may also be due to the 
fact that Zp has larger one- loop finite corrections, ~ 35%, and that perhaps 
the higher order terms, which are not accounted for by one loop boosted 

^°The plateau could eventually appear more clearly at larger values of (3. 
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perturbation theory, are important. 

We also present, in fig. H, Zs/Zp as a function of iJ?. Tlie ratio Zs/Zp 
lias a befiaviour similar to Z^l. Some discretization effects appear to be 
smaller than those observed in Zs and Zp separately. The numerical results 
support, qualitatively, the theoretical interpretation. In the intermediate 
range of /x, this ratio stays almost constant. At small or large values of 
yU the results are unstable, due to non-perturbative or discretization effects 
respectively. At ^^o? ~ 1, the value of this ratio obtained from our non- 
perturbative method is in reasonable agreement with that determined with 
the Ward identity method, Zg/Zp = 1.64(5) §Q 



7 Conclusion 

We have proposed a new non-perturbative method to renormalize lattice 
composite operators. It can be used in all cases and is particularly useful 
when the Ward identity method is not applicable. This method avoids the 
need to perform any calculations using lattice perturbation theory in the 
computation of physical quantities from lattice simulations. The success of 
our proposal is subject to the existence, at current values of P, of a window 
between the non-perturbative region at low momenta and the region of large 
momenta, where discretization errors become important. This limitation is 
however, common to all methods (with the exception of the Ward identity 
method, which can only be applied to a few cases, corresponding to finite 
operators). The window, necessary to implement the renormalization pro- 
gramme described here, seems to exists already at /3 = 6.0 and we expect that 
the range of momenta, useful for non-perturbative renormalization, will be- 
come larger as (3 increases. We are planning to apply the approach presented 
in this paper to the renormalization of the AS" = 2 four fermion operator 
(|T3) and of the heavy-light axial vector current in the static theory. 
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The value published in ref. [BJ, using half of the present stastistics, was Zp/Zs 



0.64(2), corresponding to Zs/Zp = 1.56(5). 
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Appendix 

In this appendix we demonstrate that the non-perturbative contributions to 
forward two-quark matrix elements are suppressed as inverse powers of p^, 
where p is the external momentum. Let us consider the non-amputated, 
two-fermion Green function 



Foip) = J d'xd'ye-'P-^q{0)rq{x)Or{y)), 



where Or{y) = 'i{y)^l{y) and V is one of the Dirac matrices. It is natural to 
worry about the infrared behaviour of Fo{p), because it contains the inser- 
tion of a zero-momentum operator. On general grounds, we expect infrared 
divergences to appear in Green functions computed at exceptional momenta. 
In physical terms, this is due to the fact that the matrix elements of Or, at 
zero momentum, carry information about the physical mass spectrum, which 
is unaccessible to perturbation theory. In particular, if the operator Or is 
the pseudoscalar density (75) and we work in the chiral limit, a contribution 
from a Goldstone pole, corresponding to the pion, should appear in Fo{p)- 
The presence of non-perturbative contributions to Fo{p) is signaled by the 
appearence of infrared divergences in the perturbative expansion. This is 
what happens in general. For quark degrees of freedom in the limit of large 
p^ however, the Operator Product Expansion (OPE) guarantees that Fo{p) 
can be reliably computed at all orders in perturbation theory. The argu- 
ment goes as follows. In the limit of large p^, the dominant contribution to 
Fo{p) comes from regions of integration (in x and y) where the integrand is 
singular. This implies that, in this limit, we have to study the behaviour of 

Fo{x,y) = {q{0)'^<l{x)Or{y)), 

when a; ~ and x ~ y. In both cases of course, we can use the OPE in the 
following form0 

qa{x)Or{y) -^ AlJx - y)Tpsqs{y) + • • • , 



^^From a technical point of view, OPE is a weak operator statement and, when Or is 
inserted in a Green function containing other operators, a shghtly modified form would be 
appropriate. It is easy to convince oneself, however, that the additionals terms, required 
in the general case, are not relevant to the present discussion. 
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for X r^ y, and 

g(x)rg(0)^5r(x2)g(0)rg(0) + ..., 

for a; ~ 0. 

Since QCD is asymptotically free, the behaviour oi A^p{x — y) and Br{x) 
can be computed in perturbation theory. The interaction gives rise to loga- 
rithmic corrections to the free field behaviour 



2^2^ 



Ar{x) = :/:Cr{x') = ^ 



ar\n'^^{fi^x 



x^ 



5r(a;2) = 6rln*(/iV) + ... 
Using the above equations, in the limit p^ ^ oo, we find 

'ci^xdV-'^-" ((g(O)rAr(x - y)Tqiy)) + Brix'){qiO)TqiO)Oriy)) 



where 

SM =JdSe-'P-y{m'^Yrq{y)) MO) = Jd'y{m'^q{0)Or{y)). 

The term —ip ^Cy{p'^)Sy{p) is computable in perturbation theory. On dimen- 
sional ground, for p^ — » oo, it behaves as 

• n ( '2\Q^^( \ In^'^bV/^^) 

-W^.CT[P )^r(p) ^cr- 



p2 



As for Brip"^), we find 

ln^r(pV/i2 



Br{p') ^ dr- 



p4 

By{p^) is, however, multiplied by Ar(0), which is a non-perturbative quan- 
tity. We have shown however, that, in the large p^ limit, the perturbative 
contribution dominates by one power of p^ over the non-perturbative one. 
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This is the reason why, even at an exceptional external momentum, Fo{p) is 
infrared safe in perturbation theory, when p"^ is large. It is straightforward 
to trace the suppression of the infrared divergence, by analyzing the corre- 
sponing Feynman diagrams. For example, the infrared divergent part of the 
diagram in fig. ^a is suppressed by one power oi p^, due to the momentum 
flow through the lower gluon line. 

It is interesting to see what would happen in a (fictitious) theory, with 
quarks represented by scalar fields. In this case the infrared divergence is not 
supressed, as can be seen from the diagram in fig. ^-b. The general argument 
goes as follows 



F(/) = J dSdV-^^-^(0(O)0(a;) J(y)), 

with 

J{y)=^{y)<P{y). 
When a; ~ and x ^ y, the OPE gives 

0(x) j(y) -^a[{x- yf) 4>{y) + ... mm ^ B{x'')mm + 

with 

= «ln"(/iV) + ... ^ , ^ Mn^(/xV) + . . . , 

X'^ 

and we would get 

F(p2) _ A(p2)^(/) + 5(p2) A(0) 
with 

s{p') = fd'ye-'p-y{mHy)) m = fd'yimmAy))- 



Mp')S{p' 



One now obtains 



and 



ln^(pV/i2) 



pi 



B{p^) -^ d 



ln^(pV/i^ 
p4 



The latter equations imply that non-perturbative and perturbative contri- 
butions have the same p"^ dependence. Thus, we expect the appearance of 
infrared divergences in F{p'^). 
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Figure 4: ZyL as a function of fi^a"^. Since the lattice is not symmetric 
in space and time, the results obtained by using the time- component of the 
local current (Vq) and the space components, averaged in the three directions, 
(Vk) are shown separately. The dashed line is Zyf^^ from BSPT, and the 
curve is from BDPT, on a volume of size 16^ x 32, see sec. ^, see sec. [^. 
The straight (continous) line is the result obtained using the Ward identities 
method ^, ^. 
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Figure 5: Zs as a function of fi'^a'^. We give the renormalization constant 
obtained by using the two possible definitions of the quark wave function 
renormalization, denoted by Z^, eq. ^1[), and Z'^, eq. ^C\). The dashed curve 
is Zf^^""", from boosted standard perturbation theory in the infinite volume 
limit, and the full curve is the result obtained in BDPT. 
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Figure 6: Zj^l as a function of fi'^a'^. We give the renormalization constant, 
obtained by using the time component Aq and the space components Ak, sepa- 
rately. The dashed line is Z^^^^ , from boosted standard perturbation theory, 
and the full curve is from BDPT. The straight (continous) line is the result 
obtained using the Ward identities method ^ \2\] . 
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Figure 7: Zp as a function of /x^a^ . We give the renormalization constant 
obtained by using the two possible definitions of the quark wave function 
renormalization, denoted by Z^p, eq. ^1[), and Z'^, eq. ^C\). The dashed curve 
is Zp^^"^, from boosted standard perturbation theory in the infinite volume 
limit, and the full curve is the result obtained in BDPT. 



34 



Figure 8: Zs/Zp as a function of ira . The dashed line corresponds to 
BSPT, the full line comes from the determination of this ratio, using the 
Ward identity method ^. 
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Figure 9: Diagrams with standard, a), and scalar quarks, b). 
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